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I.  INTRODUCTION 

In  this  report/  we  shall  develop  constitutive  relations 
for  fluid-saturated  porous  media/  suitable  for  inclusion  in 
standard  hydrodynamic  codes  (e.g./  CRAM  or  SKIPPER).  The 
theoretical  formulation  is  based  on  the  models  for  fluid- 
saturated  rock  aggregates  previously  developed  by  Garg  and 
Nur  [1973]  and  Garg,  et  al . [1975].  In  the  present  analysis, 
we  shall  assume  that  there  is  no  relative  motion  between  the 
fluid  and  the  solid  (no  fluid  diffusion) ; this  assumption  is 
equivalent  to  requiring  that  [Garg,  et  al . , 1975] 


<P2  uf 

where 

X - length  scale  of  interest 
Uf  - fluid  viscosity 
P - density  of  porous  media 

o 

C^  - speed  of  sound  in  porous  media 
k - permeability  of  porous  media 
$ - porosity 

It  is  straightforward  to  verify  that  the  above  inequality  is 
satisfied  for  many  field  situations  involving  shock/seismic 
wave  propagation  in  fluid-saturated  porous  media.  The  theo- 
retical formulation  will  be  outlined  in  Section  II.  A 
possible  method  for  including  the  new  constitutive  model 
into  standard  hydrodynamic  codes  is  described  in  Section  III. 
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II.  THEORETICAL  FORMULATION 

"No  fluid  diffusion"  hypothesis  implies  that  the  fluid 
mass  M in  a volume  V of  the  rock/fluid  composite  remains 
constant.  Mathematically, 


Noting  that 

Vf  * volume  occupied  by  the  fluid  - $V  , (2a) 

and  defining 

Pf  * fluid  density  = M/Vf  , (2b) 

we  have  from  (1) 

4>  v 

Pf'P,f*iV1'  (3) 

where  subscript  0 denotes  the  value  of  the  subscripted 
variable  at  t « 0.  For  a partially  saturated  media,  fluid 
density  as  defined  by  (2b)  is  not  equal  to  the  real  density 
of  the  fluid.  (In  other  words,  we  replace  the  real  fluid 
by  an  extended  fluid  occupied  by  the  volume  Vf.)  This  is 
particularly  convenient  for  treating  an  initially  partially 
saturated  medium  which  upon  shock  loading  (and  resulting 
pore  compaction)  becomes  fully  saturated;  in  this  case  the 
fluid  pressure  pf  will  remain  zero  as  long  as  pf  is  less 

than  some  reference  density  (say  p^) . 

0 

We  now  introduce  the  following  mixture  (solid/fluid) 
quantities: 

P (mixture  density)  - (l-#)p  + <pp.  (4a) 

S £ 


(mixture  stress)  * -pc  + s^.  (4b) 
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Pc  (mixture  pressure)  = (l-<J>)pg  + $pf  (4c) 

(!-<())  P E + <))fpEf 

E (mixture  internal  * ^-= = — =.  , (4d) 

energy/unit  mass)  p 

where 

p = rock  grain  density 

s 

Sij  *=  mixture  deviatoric  stress  tensor 

Pg(Pf)  “ solid  (fluid)  pressure,  and 

Eg(E^)  * solid  (fluid)  internal  energy  per  unit  mass. 

The  mixture  particle  velocity  v is  identically  equal  to  the 
solid  (fluid)  velocity  v (v,). 

— S “I 

Mixture  density  p,  velocity  v,  and  internal  energy 
E are  governed  by  the  usual  balance  equations  for  a single 
continuum: 

§£  + P div  v - 0 (5a) 

Dv 

P Dt  “ div  o + pg.  (5b) 

DE 

p Dt  “ 2 1 ' (5c) 

where 

Dt  lt  + v-  7Z  * (5c) 

We  will  define  the  mixture  strain-rate  tensor  to  be 
the  symmetric  part  of  the  velocity  gradient.  Thus 

S - y t»v  + (VJt)  (6) 
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The  strain-rate  tensor  may  be  decomposed  into  isotropic  (e) 
and  deviatoric  (e)  parts  as  follows: 


t i + e 


where  I is  the  unit  tensor. 


Mixture  density  p is  related  to  e through  the  re- 


lation: 


P * P exp(-e) 
0 


Also,  we  have 


V /V 
o 


exp(-e) 


Given  e (Eq.  6)  and  <j>  (see  below)  , p,  pf  and  pg  may 
then  be  determined  from  Eqs.  (8),  (3)  and'  (4a)  respectively. 

Now,  we  shall  consider  the  question  of  energy 
partitioning  between  the  solid  and  fluid  phases.  Let  t^ 
and  tr  denote  the  time  of  interest  and  the  thermal  relaxa- 
tion time  respectively.  If  ti  <<  tr,  the  adiabatic  assump- 
tion (no  thermal  exchange  between  the  solid  and  fluid  phases) 
may  be  employed.  Similarly,  the  isothermal  assumption 
(Tf  * Tg)  is  appropriate  in  the  reverse  case,  i.e.,  when 
>>  tr.  In  the  intermediate  range  (tr  - Ott^),  it  is 
necessary  to  adopt  a constitutive  model  for  heat  flow  be- 
tween the  two  phases  (see,  e.g. , Riney,  et  aL  [1971]).  In 
laboratory  situations,  the  adiabatic  assumption  is  generally 
adequate.  Adiabatic  assumption  may  not,  however,  be  always 
appropriate  in  field  applications.  For  the  sake  of  sim- 
plicity, we  will  adopt  the  adiabatic  assumption  in  the 
present  analysis.  (It  cannot  be  over-emphasized  that  this 
assumption  may  be  invalid  in  specific  applications.)  In 
this  case,  the  fluid  internal  energy  can  be  written  as: 


* 


R-2766 


t 

Ef  ° • / pf d (57)  (10> 

J0  t 

Given  Ef  (Eq.  10)  and  E (Eq.  5c),  Eg  may  be  deter- 
mined from  Eq.  (4d)  . 

We  now  wish  to  prescribe  the  constitutive  relations 
for  the  pressures  (pg,  pf)  and  deviatoric  stresses  S.  In 
the  rock— fluid  mixture,  only  the  rock  can  sustain  shear 
stresses;  and  one  can  therefore  directly  postulate  shear 
laws  for  the  entire  composite.  We  will  assume  the  solid  to 
be  elastic-plastic;  the  shear  stresses  are  assumed  to  he 
restricted  by  a yield  criterion.  The  shear  laws  are,  there- 
fore,  as  follows; 

(i)  Elastic  response  or  for  unloading  from  yield 
surface  (S  : S < 2Y2): 


(ii) 


2ue 


(11) 


Here  \i  denotes  the  shear  modulus  of  the 
porous  rock  and  Y is  the  yield  stress  in 
shear. 


Plastic  response; 
S ; S - 2Y2 


(12) 


where 
Y = Y(P) 


(13a) 


1/3 


Pc  " Pf 


k(±)  . 


and 


(13b) 


* third  invariant  of  deviatoric  stresses 
- Det  (S)  . 


(13c) 


f 
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In  writing  down  Eq.  (13a) , we  have  assumed  that  the  yield 
surface  Y of  a porous  rock  is  governed  by  the  "effective 
stress  law"  (see,  e.g.,  Garg  and  Nur  [1973]).  A mixture 
model  for  pressures  (pg,  pf)  can  be  formulated  by  assuming 
that  the  pressure  law  (p  versus  p and  E relationship) 
for  each  constituent  as  a single  continuum  applies  in  the 
mixture.  Thus 

Pf  = Pf  (Pf  / Ef ) 


= PS(P3' 


Es> 


(14b) 


The  above  constitutive  formulation  is  complete  only 
when  porosity  $ is  prescribed.  Porosity  $ will,  in 
general,  depend  upon  pc,  p^,  S,  Eg,  E^,  and  the  past 
loading  history.  In  shock  wave  propagation  th  rough  dry 
rocks,  it  is  often  assumed  (perhaps  without  much  justifica- 
tion) that  4>  depends  only  upon  pressure  pc-  For  wet 
rocks,  a similar  model  is  formulated  by  assuming  that  <$> 
is  a function  of  only  pc  and  p^  (see  Gc.rg,  et  al , [1975]). 

0 - 4>(PC*  Pf)  (1 


Empirical  observations  under  isotropic  loading  reveal  that 
<p  is  a particularly  simple  function  of  pc  and  p^;  <P 
depends  to  a close  approxi  -.ation  on  Pc  “ Pf  • 

<♦>  * 4>(PC  “ Pf) 


Function  <p  (pc  - pf ) may  be  determined  from  hydrostatic 
tests  with  pf  kept  at  zero.  This  completes  the  theoretical 
formulation.  In  the  next  section,  we  shall  outline  an 
iterative  procedure  for  incorporating  this  model  into  conven- 
tional hydrodynamic  codes. 
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III.  ITERATIVE  PROCEDURE  FOR  INCORPORATING 
THE  THEORETICAL  MODEL  INTO 
COMPUTER  CODES 


In  this  section,  we  will  describe  a possible  itera- 
tive procedure  for  incorporating  the  theoretical  model  of 
Section  II  into  standard  hydrodynamic  codes.  This  pro- 
cedure has  not  oeen  tested,  and  may  require  minor  adjust- 
ments to  ensure  convergence.  For  the  sake  of  convenience, 
the  iterative  procedure  is  presented  below  in  a flow  chart 
form.  (In  the  following,  superscript  n denotes  the  value 
of  the  superscripted  variable  at  time  t = nAt) . Given 
strain-rate  tensor  c at  time  t * (n+1)  At  and  Pp,  <f>n. 


,n 


Pc- 


*f' 
solves  for 
„n+l 


_n 
P / 


Sn 
2 • 


n+1 
Pf  , 


en 

.n+1 

# 


and 


E", 


the  iteration  procedure 

;.n+l  n+1  n-*-l  n+1  . 

t P„  , P , S , and 


Iteration  Procedure 


(Saved  arrays  are  pf,  <p,  Ef,  pc,  s,  E) 


Step  1.  Initialize  variables 


n+1 


n 


_n+l  n 
Pf  » Pf 


.n+1  .n 
<p  « 4> 


E^+1  - (pnEn  - pJ<t>nE5]/[p£(l-$n)] 
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Stop  2. 


Determine 
e ■ Trace 


e ■ e 

A*  a» 

_n+l 


n+1 


new  strains 
e 

a* 

I 

A# 

+ At  £ 

+ At  e 

A# 


Step  3. 
Step  3A. 


Start  of  Iteration  Loop 


n+1  n+1  .n+1 


» » 


Calculate  p 

pn+1  = p exp(-£n+1) 

0 

p T 1 - °.£  “p(-^n+1> 


s 


»fl  - tp«+i . *"+V£+1>/<iV+1> 


Step  3B.  Calculate  new  shear  stresses 
w - Max  ( [ p (1-E*J+1/E  )J,  0) 

0 S Ml 

- shear  modulus  at  room  temperature) 

(Em  - specific  internal  energy  for  the  solid 
phase  at  incipient  melting) 

AS  ■ 2y  At  e 

A*  A# 

Sn+1  « sn  + AS 

A*  A*  A# 

J'n+1  - [Sn+1  s Sn+1]/2 
2 - * 

J'n+1  . Det [Sn+1] 

i ~ 
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pH+1 


/ lT  \ 

z^n+1  _n+l,  . 1 / J3  \ 

(pc  " Pf  ) + l(“V-/ 


1/3 


Yn+1  = Y (Pft+1)  x (1  - En+1/E  ) 
o s m 


Yn+1  * Max (Yn+1 , 0) 


( Y ^ (P ) - yield  strength  at  room  temperature) 


If  Vj-+1  > Yn+1,  then  set  Sn+1  - -n+1  Y 


n+1 


H+T 


Step  3C.  Calculate  new  internal  energies 


v = -p/p2  « e/p 


n+1 


E « -(p  + q)n+1v  + Sn+1  : e 


(q  - artificial  viscosity 
v - specific  volume  » 1/p) 

En+1  = En  + At  E 


Avf  - A 


kn+l  .n 


(v^  - specific  vojume  of  fluid) 


AEf  « -Pf+1  AVj 


Ef+1  " Ef  + AEf 


,n->  1 


ipn+1En+1  - efV+1E£+1i/cP;;+1<i-*n+1>] 
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IV.  CONCLUDING  REMARKS 


The  preceding  sections  describe  a simple  constitu- 
tive model  for  fluid-saturated  porous  rocks.  In  developing 
the  constitutive  model,  it  was  assumed  that  (1)  there  is  no 
relative  motion  between  the  fluid  and  the  solid,  (2)  no  sig- 
nificant heat  exchange  occurs  between  the  solid  and  fluid 
phases,  and  (3)  porosity  ♦ is  a function  of  only  the  mix- 
ture pressure  pc  and  fluid  pressure  pf.  The  last  two 
assumptions  may  not  be  valid  under  certain  field  conditions. 
These  assumptions  can,  however,  be  relaxed  at  the  expense  of 
some  simplicity  [Riney,  et  ah,  1971;  Garg,  1975]. 
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